
library(readstata13)
library(ggplot2)
library(stats)
library(dplyr)

rm(list=ls())
setwd("replication_APSR")

#######multiplot function from -http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

#######################################################
#######################################################
# Councilor perfromance composite index
#######################################################
#######################################################

##################################
# Kernel distribution of councilor perfromance composite index
# Figure 21 online appendix
##################################
MainOutcomes=read.dta13("Analysis_final/MainOutcomes.dta", convert.factors = TRUE, nonint.factors=TRUE)
names(MainOutcomes)
head(MainOutcomes)

(x=mean(MainOutcomes$index_PerformanceA, na.rm=TRUE))
pdf(file="figures/CompositeIndexKernelDist.pdf", width=12, height=6) 
ggplot(MainOutcomes, aes(index_PerformanceA)) + geom_density(colour="darkgreen", size=2, fill="green")+ ggtitle("Composite Index Kernel Distribution") + theme_bw(base_size = 12, base_family = "Helvetica")+
  scale_x_continuous(name="Councilor performance composite index")+geom_vline(xintercept = x, colour="red") + 
  theme(panel.grid.major.y = element_line(),panel.grid.major = element_line(), 
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(size = 18, face = "bold", vjust = 1.5),axis.title.y = element_text(colour="grey20",size=14,face="bold"),axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"),  axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank())
dev.off()

##########################################################################
# Box plot of perfromance composite index by treatment and electoral pressure 
# Figure 22
##########################################################################
IndPerBox= MainOutcomes %>%  filter(!is.na(ID), !is.na(competitive), !is.na(index_PerformanceA)) %>%select(id, ID,competitive, index_PerformanceA)

pdf("figures/CompositeIndexTreat.pdf", width=12, height=7)
ggplot(IndPerBox, aes(y=index_PerformanceA, x=ID, fill=ID)) + geom_boxplot() + scale_y_continuous(name="Complosite Index") + 
  geom_jitter(alpha=0.5, aes(color=ID),position = position_jitter(width = .2)) +  coord_flip() + facet_grid(competitive ~ .)+
  theme_bw() + scale_x_discrete(name="") + scale_fill_grey(start=.5, end=1) + stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=4,show.legend = FALSE) +
  theme(plot.title = element_text(size = 20, face = "bold", vjust = 1.5), axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=14,face="bold"), axis.text.y = element_text(colour="grey20",size=14,face="bold"), legend.position="none",  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()) +ggtitle("Composite Index by treatment")
dev.off()

#######################################################
#######################################################
# ACODE scorecard
#######################################################
#######################################################

#######################################################
# ACODE baseline scorecard by councilor mandate (violin)
#Figure 10 online appendix
#######################################################
demo=read.dta13("Analysis_final/Councilors_covariates_v2.dta", convert.factors = TRUE, nonint.factors=TRUE)
BaseScore= demo %>%  filter(!is.na(ID), !is.na(competitive), !is.na(total_score_2011_2012)) %>%select(id, ID,competitive, total_score_2011_2012, mandate)
head(BaseScore)

range(BaseScore$total_score_2011_2012, na.rm=T)

pdf(file="figures/Score2011Violin.pdf", width=12, height=6) 
ggplot(BaseScore, aes(factor(mandate), total_score_2011_2012)) + geom_violin(aes(fill = mandate)) + geom_boxplot(width=.1)+
  ggtitle("ACODE's Total Score at Baseline") + theme_bw(base_size = 12, base_family = "Helvetica")+
  scale_y_continuous(name="Total Score 2011-2-12")+ scale_x_discrete(name="") + coord_flip() +
  theme(panel.grid.major.y = element_line(),panel.grid.major = element_line(), legend.position = "none",
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(size = 18, face = "bold", vjust = 1.5),axis.title.y = element_text(colour="grey20",size=14,face="bold"),axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"),  axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank())
dev.off()

##################################
# Distribution of last total score (2014-2015)
#Figure 11 online appendix
##################################
MainOutcomes=read.dta13("Analysis_final/MainOutcomes.dta", convert.factors = TRUE, nonint.factors=TRUE)
head(MainOutcomes)
(x=mean(MainOutcomes$total_score_2014_2015, na.rm=TRUE))

pdf(file="figures/TotaScoreDensity.pdf", width=12, height=6) 
ggplot(MainOutcomes, aes(total_score_2014_2015)) + geom_density(colour="darkgreen", size=2, fill="green")+ ggtitle("Total Score Distribution") + theme_bw(base_size = 12, base_family = "Helvetica")+
  scale_x_continuous(name="ACODE Total Score (2014-2015)")+geom_vline(xintercept = x, colour="red") + 
  theme(panel.grid.major.y = element_line(),panel.grid.major = element_line(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(size = 18, face = "bold", vjust = 1.5),axis.title.y = element_text(colour="grey20",size=14,face="bold"),axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"),  axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank())
dev.off()

######################################################
# Boxplot of 2013-2014 score by treatment 
#Figure 13  online appendix
######################################################
ACODEScores= MainOutcomes %>% select(id, ID,competitive, total_score_2014_2015) %>%  
  filter(!is.na(ID), !is.na(competitive),!is.na(total_score_2014_2015))

pdf("figures/ACODEScoresTreat.pdf", width=12, height=7)
ggplot(ACODEScores, aes(y=total_score_2014_2015, x=ID, fill=ID)) + geom_boxplot() + scale_y_continuous(name="Peers' evaluation") + 
  geom_jitter(alpha=0.5, aes(color=ID),position = position_jitter(width = .2)) +  coord_flip() + facet_grid(competitive ~ .)+
  theme_bw() + scale_x_discrete(name="") + scale_fill_grey(start=.5, end=1) + stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=4,show.legend = FALSE) +
  theme(plot.title = element_text(size = 20, face = "bold", vjust = 1.5), axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=14,face="bold"), axis.text.y = element_text(colour="grey20",size=14,face="bold"), legend.position="none",  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()) +ggtitle("Peers' evaluation by treatment")
dev.off() 

##################
#Figure 9 online appendix
##################

scores=read.dta13("estimates/ScoresByTreatment.dta", convert.factors = TRUE, nonint.factors=TRUE)
head(scores)
colnames(scores)[1] <- "Treat"

pdf("figures/TotalScoreOverTime.pdf", width=12, height=7)
ggplot(scores, aes(x=factor(year), y=total, fill=Treat, group=Treat, linetype=Treat)) + 
  geom_line(size=1.5) + geom_point(size=4.5, shape=21) + facet_grid(.~competitive)+
  scale_x_discrete(name="") + theme_bw() + ggtitle("ACODE's Total Score by Treatment") +
  scale_y_continuous(name="Total Score", breaks = round(seq(min(40), max(60), by =2),0), limits = c(40, 60)) +
  theme(plot.title = element_text(size = 18, face = "bold", vjust = 1.5),
        axis.title.y = element_text(colour="grey20",size=15,face="bold"),
        axis.text.x = element_text(colour="grey20",size=15,face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"), 
        legend.position=c(.5, .2), legend.background = element_rect(fill="gray90", size=.5, linetype="dotted"),
        legend.key.size = unit(1.6, "cm"),
        axis.title.x = element_text(colour="grey20",size=15,face="bold"), axis.ticks.y = element_blank())
dev.off()

#######################################################
#######################################################
# Technocrats' assessment
#######################################################
#######################################################

###############################################
# Distribution of technocrat mean assessment
#Figure 14 online appendix
###############################################
TechP=read.dta13("Analysis_final/Technocrats_assessment.dta", convert.factors = TRUE, nonint.factors=TRUE)
head(TechP)
TechPM= TechP%>%  filter(!is.na(index_technocratsA))
(x=mean(TechPM$index_technocratsA, na.rm=TRUE))

Demos=read.dta13("Analysis_final/Councilors_covariates_v2.dta", convert.factors = TRUE, nonint.factors=TRUE)
TechDem=merge(TechPM, Demos, by="id")
head(TechDem)

pdf(file="figures/TechIndexKernelDist.pdf", width=12, height=6) 
ggplot(TechPM, aes(index_technocratsA)) + geom_density(colour="darkgreen", size=2, fill="green")+ ggtitle("Distribution of technocrats' mean assessement") + theme_bw(base_size = 12, base_family = "Helvetica")+
  scale_x_continuous(name="technocrats' mean assessement (std)")+geom_vline(xintercept = x, colour="red") + 
  theme(panel.grid.major.y = element_line(),panel.grid.major = element_line(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(size = 18, face = "bold", vjust = 1.5),axis.title.y = element_text(colour="grey20",size=14,face="bold"),axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"),  axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank())
dev.off()

######################################################
# Technocrats' eveluation by treatment 
# Figure 15 online appendix
######################################################

pdf("figures/TechMeanEvalTreat.pdf", width=12, height=7)
ggplot(TechDem, aes(y=index_technocratsA, x=ID, fill=ID)) + geom_boxplot() +  scale_y_continuous(name="Technocrats' assessement") + 
  geom_jitter(alpha=0.5, aes(color=ID),position = position_jitter(width = .2)) + coord_flip() + facet_grid(competitive ~ .)+
  theme_bw() + scale_x_discrete(name="") + scale_fill_grey(start=.5, end=1) + stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=4,show.legend = FALSE) +
  theme(plot.title = element_text(size = 20, face = "bold", vjust = 1.5), axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=14,face="bold"), axis.text.y = element_text(colour="grey20",size=14,face="bold"), legend.position="none",  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()) +ggtitle("Technocrats' assessement summary index by treatment")
dev.off() 

################################################
# Technocrats mean assessment correlation with ACODE's total scorecard score
# Figure 20 online appendix
################################################
names(MainOutcomes)
(techcor=round(cor(MainOutcomes$total_score_2013_2014, MainOutcomes$technocrats, use = "complete"),3))

pdf("figures/TotalScoreMeanTech.pdf", width=12, height=7)
(tech=ggplot(MainOutcomes, aes(x=technocrats, y=total_score_2013_2014)) + geom_point(size=3) +
  geom_smooth(method=lm, color="darkred", fill="blue") + theme_bw() +  ggtitle("Technocrats' mean assessement vs. Total score") + 
  scale_y_continuous(name="Total Score (2013-2014)", limits=c(0, 80))  + scale_x_continuous(name="Technocrats' mean assessement")+
  geom_text(aes(x=2, y=64, label=paste("r^2==", techcor), group=NULL), size = 6, color = "red", data=MainOutcomes, parse = T) +
  theme(plot.title = element_text(size = 18, face = "bold", vjust = 1.5),
        axis.title.y = element_text(colour="grey20",size=15,face="bold"),
        axis.text.x = element_text(colour="grey20",size=13,face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"), 
        axis.title.x = element_text(colour="grey20",size=15,face="bold"), axis.ticks.y = element_blank()))
dev.off() 

fit <- lm(total_score_2013_2014 ~ technocrats, data=MainOutcomes)
summary(fit)

#######################################################
#######################################################
# Peers' mean evaluation
#######################################################
#######################################################

#######################################################
# Peers' mean evaluation: Kernel density Distribution
# Figure 12 online appendix
######################################################

names(MainOutcomes)
(x=mean(MainOutcomes$avescore, na.rm=TRUE))

pdf(file="figures/AveScoreStdKernelDist.pdf", width=12, height=6) 
(kernel=ggplot(MainOutcomes, aes(avescore)) + geom_density(colour="darkgreen", size=2, fill="green")+ ggtitle("Peer mean evaluation kernel distribution") + theme_bw(base_size = 12, base_family = "Helvetica")+
  scale_x_continuous(name="Peer mean evaluation")+geom_vline(xintercept = x, colour="red") + 
  theme(panel.grid.major.y = element_line(),panel.grid.major = element_line(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(size = 18, face = "bold", vjust = 1.5),axis.title.y = element_text(colour="grey20",size=14,face="bold"),axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"),  axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()))
dev.off()

######################################################
# Boxplot of mean peers' eveluation by treatment 
######################################################
PeerScores= MainOutcomes %>% select(id, ID,competitive, avescore, total_score_2013_2014) %>%  
  filter(!is.na(ID), !is.na(competitive),!is.na(total_score_2013_2014), !is.na(avescore))

pdf("figures/PeersMeanEvalTreat.pdf", width=12, height=7)
ggplot(PeerScores, aes(y=avescore, x=ID, fill=ID)) + geom_boxplot() + scale_y_continuous(name="Peers' evaluation") + 
  geom_jitter(alpha=0.5, aes(color=ID),position = position_jitter(width = .2)) +  coord_flip() + facet_grid(competitive ~ .)+
  theme_bw() + scale_x_discrete(name="") + scale_fill_grey(start=.5, end=1) + stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=4,show.legend = FALSE) +
  theme(plot.title = element_text(size = 20, face = "bold", vjust = 1.5), axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=14,face="bold"), axis.text.y = element_text(colour="grey20",size=14,face="bold"), legend.position="none",  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()) +ggtitle("Peers' evaluation by treatment")
dev.off() 

####################################################################
# Peer evaluation: Box-plot of total score (using data in long format)
# Figure 19 (right) online appendix
#####################################################################
PeerL=read.dta13("Analysis_final/Peer_evaluations_Long.dta", convert.factors = TRUE, nonint.factors=TRUE)
names(PeerL)

PeerLDem=merge(PeerL, PeerScores, by="id")
PeerLDem$pscore=as.numeric(PeerLDem$peer_score)
names(PeerLDem)
dim(PeerLDem)

peerslongM= PeerL%>%  filter(!is.na(total_score_2013_2014), !is.na(pscore))
head(peerslongM)
dim(peerslongM)

(peerslongF=ggplot(peerslongM, aes(y=total_score_2013_2014, x=peer_score, fill=peer_score)) + geom_boxplot() + scale_y_continuous(name="Total Score (2013-2014)") + 
  geom_jitter(alpha=0.5, aes(color=peer_score),position = position_jitter(width = .2)) +
  theme_bw() + scale_x_discrete(name="") + scale_fill_grey(start=.5, end=1) + stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=4,show.legend = FALSE) +
  theme(plot.title = element_text(size = 20, face = "bold", vjust = 1.5), axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=14,face="bold", angle = 30, hjust = 1), axis.text.y = element_text(colour="grey20",size=14,face="bold"), legend.position="none",  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()) +ggtitle("Total Scores | Peer Assessement"))

#####################################################
# peer mean evaluation correlation with total score
# Figure 19 (left) online appendix
#####################################################
peerslongM2= peerslongM%>%  filter(!is.na(total_score_2013_2014), !is.na(avescore.x))

peerswide=peerslongM2[peerslongM2$tag==1,]
names(peerswide)
peerswide$avescore=as.numeric(peerswide$avescore.x)
peercor=round(cor(peerswide$total_score_2013_2014, peerswide$avescore, use = "complete"),3)

(peer=ggplot(peerswide, aes(x=avescore, y=total_score_2013_2014)) + geom_point(size=3) + geom_smooth(method=lm, color="darkred", fill="blue") +
  scale_x_continuous(name="Peer (mean) assessement") + theme_bw() +
  ggtitle("Total Score | Peer Mean Assessement") + scale_y_continuous(name="Total Score (2013-2014)",limits=c(0, 80)) +
  geom_text(aes(x=1.5, y=71, label=paste("r^2==", peercor), group=NULL), size = 6, color = "red", data=peerswide, parse = T) +
  theme(plot.title = element_text(size = 18, face = "bold", vjust = 1.5),
        axis.title.y = element_text(colour="grey20",size=15,face="bold"),
        axis.text.x = element_text(colour="grey20",size=13,face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"), 
        axis.title.x = element_text(colour="grey20",size=15,face="bold"), axis.ticks.y = element_blank()))

pdf("figures/TotalScoreMeanPeer.pdf", width=12, height=7)
multiplot(peer, peerslongF, cols=2)
dev.off() 

#######################################################
#######################################################
# School grant application
#######################################################
#######################################################

#######################################################
#  Kernel density Distribution of complete school grants
# Figure 17 online appendix
#######################################################
SCBox= MainOutcomes %>%  filter(!is.na(ID), !is.na(competitive), !is.na(CompleteApp), !is.na(AllApp)) %>%select(id, ID,competitive, AllApp,CompleteApp)
(x=mean(SCBox$CompleteApp, na.rm=TRUE))

pdf(file="figures/AppsSchoolCompleteKernelDist.pdf", width=12, height=6) 
(kernel=ggplot(SCBox, aes(CompleteApp)) + geom_density(colour="darkgreen", size=2, fill="green")+ ggtitle("School Grant Applications Distribution") + theme_bw(base_size = 12, base_family = "Helvetica")+
  scale_x_continuous(name="Complete School Grant Applications")+geom_vline(xintercept = x, colour="red") + 
  theme(panel.grid.major.y = element_line(),panel.grid.major = element_line(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        plot.title = element_text(size = 18, face = "bold", vjust = 1.5),axis.title.y = element_text(colour="grey20",size=14,face="bold"),axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=13,face="bold"),  axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()))
dev.off()  

######################################################
# Box plot of complete school grants by treatment 
# Figure 18 online appendix
######################################################

pdf("figures/AppsSchoolCompleteTreat.pdf", width=12, height=7)
ggplot(SCBox, aes(y=CompleteApp, x=ID, fill=ID)) + geom_boxplot() + scale_y_continuous(name="Complete School Grant Applications") + 
  geom_jitter(alpha=0.5, aes(color=ID),position = position_jitter(width = .2)) +  coord_flip() + facet_grid(competitive ~ .)+
  theme_bw() + scale_x_discrete(name="") + scale_fill_grey(start=.5, end=1) + stat_summary(fun.y=mean, colour="red", geom="point", shape=18, size=4,show.legend = FALSE) +
  theme(plot.title = element_text(size = 20, face = "bold", vjust = 1.5), axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=14,face="bold"), axis.text.y = element_text(colour="grey20",size=14,face="bold"), legend.position="none",  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), axis.ticks.y = element_blank()) +ggtitle("Complete School Grant Applications by treatment")
dev.off() 




